Drought and other climate change-associated weather events could affect ecosystem services provided by soil microbial communities. Around 40% of the UK’s land cover is grassland, so it is important to understand how water stress affects microbial communities. This may help ecologists to identify strategies for promoting grassland resilience. The purpose of this project was to analyse 16S and ITS amplicon sequences from control and sheltered soil samples, collected by RestREco, to determine how short-term water stress affects the diversity, taxonomy and functions of soil bacteria and fungi. We used a well-established pipeline involving bioinformatics tools QIIME2, DADA2, and PICRUSt2 (see workflow diagram below), alongside R packages phyloseq, vegan, ggpicrust2 and FUNGuildR.
We also found that site differences masked the full impact of drought, but the above findings indicate there were measurable effects on our bacterial and fungal communities, although the response from fungi was much more muted
This page focuses on alpha and beta diversity and taxonomy.
For differential abundance results, visit: https://kerrycranfield.github.io/RestREcoDrought/restreco_diff_abund.html
For functional analysis, visit: https://kerrycranfield.github.io/RestREcoDrought/restreco_func.html
RestREco (Restoring Resilient Ecosystems) is a NERC-funded research project that aims to identify ecological restoration practices that can build resilience in ecosystems, particularly woodlands and grasslands. The initiative brings together researchers from:
RestREco studies a network of 133 ecological restoration sites across England and Scotland. The project splits its work into seven work packages, each focusing on different factors and relating them to ecosystem characteristics that arise from component interactions. It aims to identify key drivers of ecosystem development, including:
The goal is to understand how these factors influence ecosystem complexity, function, and resilience to future pressures (RestREco, 2024).
This study used targeted amplicon sequences (16S, ITS), retrieved from soil samples, to determine how drought affects soil microbial communities, specifically bacteria and fungi, in UK grasslands.
This study specifically focuses on:
We focused on the following four hypotheses. That drought:
A total of 36 soil samples were collected across six sites in England for each marker (6 per site: 3 control, 3 shelter).
After sieving and extracting DNA, the V4-V5 region of the bacterial 16S gene and the ITS1-1F region of the fungal ITS gene was sequenced by Novogene.
Some fungal samples failed the initial quality control during the sequencing process, resulting in the removal of one site from the analysis (for fungi only), while other sites where one or two samples failed quality checks, were kept in (also for fungi).
Sample Zone - Based on GPS coordinates
| Metric | 16S | ITS |
|---|---|---|
| Microbial group | Bacteria | Fungi |
| Region sampled | England | England |
| Number of sites | 6 | 5 |
| Samples per site | 6 | 6 |
| Total samples removed | 0 | 3 |
| Total samples | 36 | 27 |
| Total samples per treatment | 18 | 13 (for most sites) |
| Average reads per sample | ~66 000 | ~61 000 |
Sampling Summary
Before diversity, taxonomy and functional analysis could take place, the initial analysis of the DNA sequences were performed using Cranfield University’s high performance computing platform, Crescent2.
For all the scripts from the high performance computing/QIIME2 steps visit the project’s GitHub page at GitHub - RestREco Drought
These stages included quality control, trimming, denoising (merging paired reads, chimera removal and generating feature tables)
To make sure the rain shelters put in place for the experiment actually led to lower soil moisture when compared to the control, we carried out a paired Wilcoxon signed-rank test to check differences were significant.
Here, we found soil in the shelter plots was significantly drier than in the control plots.
The differences in individual sites varied, with some sites showing a much larger difference than others
Across all samples, median soil moisture content was 39.8% for control and 6% for shelter with soil moisture percentages ranging from 17.6% to 48.5% (control) and 1.5% to 34.8% (shelter)
The feature table, taxonomy table, tree and metadata were imported into phyloseq from QIIME2 using qiime2R. Sequences belonging to archaea and unassigned sequences were removed.
Below, we show how the phyloseq object changes with each step during preprocessing. A summary table is displayed after rarefaction is applied.
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 28079 taxa and 36 samples ]
## sample_data() Sample Data: [ 36 samples by 4 sample variables ]
## tax_table() Taxonomy Table: [ 28079 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 28079 tips and 28055 internal nodes ]
##
## d__Archaea d__Bacteria Unassigned
## 19 28055 5
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 28055 taxa and 36 samples ]
## sample_data() Sample Data: [ 36 samples by 4 sample variables ]
## tax_table() Taxonomy Table: [ 28055 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 28055 tips and 28031 internal nodes ]
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 6327 taxa and 25 samples ]
## sample_data() Sample Data: [ 25 samples by 4 sample variables ]
## tax_table() Taxonomy Table: [ 6327 taxa by 7 taxonomic ranks ]
##
## Fungi
## 6327
Both the bacterial and fungal feature tables were then agglomerated to genera level. Reads below this level were retained.
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1118 taxa and 36 samples ]
## sample_data() Sample Data: [ 36 samples by 4 sample variables ]
## tax_table() Taxonomy Table: [ 1118 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1118 tips and 1117 internal nodes ]
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 711 taxa and 25 samples ]
## sample_data() Sample Data: [ 25 samples by 4 sample variables ]
## tax_table() Taxonomy Table: [ 711 taxa by 7 taxonomic ranks ]
Feature tables were then aggregated to site level using speedyseq as phyloseq merge_samples combines abundances using sum. Median was felt to be more representative of sample distribution plus prevented outliers from skewing the data. Read counts were again generated to check the aggregation step was performed correctly and metadata was checked to ensure no NAs were coerced into the data.
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1118 taxa and 12 samples ]:
## sample_data() Sample Data: [ 12 samples by 4 sample variables ]:
## tax_table() Taxonomy Table: [ 1118 taxa by 7 taxonomic ranks ]:
## phy_tree() Phylogenetic Tree: [ 1118 tips and 1117 internal nodes ]:
## taxa are rows
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 711 taxa and 10 samples ]:
## sample_data() Sample Data: [ 10 samples by 4 sample variables ]:
## tax_table() Taxonomy Table: [ 711 taxa by 7 taxonomic ranks ]:
## taxa are rows
Rarefaction curves were then generated to find sampling depths using the rarecurve function in the vegan package to generate a tibble friendly data frame to be passed to ggplot2 to generate attractive rarefaction plots. As non-aggregated data would be needed for later steps (ie. alpha diversity between sites), rarefaction curves were generated for this alongside those for aggregated data.
Bacterial samples for aggregated counts were rarefied to a depth of 16 398 in order to retain all samples, while non-aggregated counts were rarefied to a depth of 13 320, using the rarefy_even_depth function with the replace parameter set to false.
For fungi, samples for aggregated counts were rarefied to a depth of 36 120 while non-aggregated counts were rarefied to 30 659.
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 627 taxa and 12 samples ]:
## sample_data() Sample Data: [ 12 samples by 4 sample variables ]:
## tax_table() Taxonomy Table: [ 627 taxa by 7 taxonomic ranks ]:
## phy_tree() Phylogenetic Tree: [ 627 tips and 626 internal nodes ]:
## taxa are rows
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1084 taxa and 36 samples ]:
## sample_data() Sample Data: [ 36 samples by 4 sample variables ]:
## tax_table() Taxonomy Table: [ 1084 taxa by 7 taxonomic ranks ]:
## phy_tree() Phylogenetic Tree: [ 1084 tips and 1083 internal nodes ]:
## taxa are rows
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 509 taxa and 10 samples ]:
## sample_data() Sample Data: [ 10 samples by 4 sample variables ]:
## tax_table() Taxonomy Table: [ 509 taxa by 7 taxonomic ranks ]:
## taxa are rows
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 708 taxa and 25 samples ]:
## sample_data() Sample Data: [ 25 samples by 4 sample variables ]:
## tax_table() Taxonomy Table: [ 708 taxa by 7 taxonomic ranks ]:
## taxa are rows
Throughout preprocessing, microbiome’s summarize_phyloseq function was used to generate statistics to see how the number of sequences changed throughout the process. These statistics included read counts and sparsity metrics (proportion of cells where count = 0).
| Metric | Start | Agglomeration | Aggregation | Rarefaction |
|---|---|---|---|---|
| Min reads | 13 222 | 13 220 | 16 398 | 16 398 |
| Max reads | 32 667 | 32 660 | 28 026 | 16 398 |
| Total reads | 865 814 | 865 546 | 276 014 | 196 776 |
| Average reads | ~24 050 | ~24 043 | ~23 001 | 16 398 |
| Sparsity | 94.6% | 66.7% | 69.2% | 45.7% |
| No. singletons | 951 | 3 | 484 | 25 |
| % singletons | 3.39 | 0.27 | 0.09 | 3.99 |
| Metric | Start | Agglomeration | Aggregation | Rarefaction |
|---|---|---|---|---|
| Min reads | 30 659 | 30 659 | 36 120 | 36 120 |
| Max reads | 78 553 | 78 553 | 56 936 | 36 120 |
| Total reads | 1 329 016 | 1 329 016 | 481 530 | 361 200 |
| Average reads | ~53 161 | ~53 161 | 48 153 | 36 120 |
| Sparsity | 91% | 91% | 70% | 59% |
| No. singletons | 0 | 0 | 203 | 10 |
| % singletons | 0 | 0 | 0.28 | 1.96 |
Bar plots showing taxonomic composition at a phylum, class and family were generated for each variable being assessed: treatment, site and group (site-treatment). Functions for each type of bar plot were created before being used to generate the plots below.
No visible shifts in the most abundant bacterial phyla were observed in response to moisture deficits. Bacterial communities were dominated by Actinomycetota, Pseudomonadota, Acidobacteriota, Bacillota I and Chloroflexota.
These phyla are common to UK grassland soils so their presence here is to be expected.
However, the effects of moisture deficits on taxonomic composition were visible within individual sites. Shifts could be seen in the relative abundances of Actinomycetota, Pseudomonadota and Bacillota I, with the direction of change being site-dependent. Actinomycetota increased while Pseudomonadota showed a slight decrease in most sites following the soil moisture reduction. One site, Harley Farms, showed the opposite trend for both taxa. For Bacillota I, changes in relative abundance between sites could be seen.
Both Actinomycetota and Bacillota I are drought-tolerant so their changes in drought conditions could be indicative of these strategies.
As at phylum level, no visible changes in taxonomic composition could be seen between treatments overall.
On a site level, no large differences in taxonomic composition could be seen. Bacilli-A (orange) appeared to vary between sites with the biggest difference between Jemma 6 and Lardon Chase.
Between treatments on each site, differences were subtle. Harley Farms showed some substantial differences between treatments for Bacilli-A (increased in drought) while Lardon Chase saw a decrease in Alphaproteobacteria in drought.
Differences in taxonomic composition between control and drought were not observed. At a site-treatment level, subtle shifts could be seen in some families (DSM-18226_301387, Geminicoccaceae for example)
Fungal community taxonomy also showed no major shifts in the abundance of the major phyla in response to water stress. Ascomycota, Mortierellomycota and Basidiomycota dominated communities. Zoopagomycota increased slightly in drought. Around 10% of fungal ASVs could not be identified to phyla level.
On individual sites, treatment-associated shifts in fungal community composition could be observed in Ascomycota and Mortierellomycota with direction of change site-dependent.
Ascomycota abundances increased in shelter samples in three out of the five sites while Mortierellomycota showed the opposite trend.
Interestingly, the sites where Ascomycota decreased under water stress were the same sites where Mortierellomycota increased. An increase in the abundance of Zoopagomycota was also of note on Windmill Farm
At a class level, no major taxonomic shifts could be seen in relation to drought. Sordariomycetes and Mortierellomycetes were the dominant classes here.
On some sites, changes in the abundances of these families could be seen between control and drought. The direction of change varied between sites, perhaps indicating an interaction between water stress and site conditions
At a family level the dominant fungal families were Mortierellaceae, Bionectriaceae and Plectosphaerellaceae. Mortierellaceae showed different directions of change in different sites.
Alpha diversity calculations were run for both bacteria and fungi to determine if there were differences in diversity between sites and treatments.
The Chao1, Shannon and Simpson metrics were used to measure richness and evenness. Chao1 is weighted towards doubletons and singletons while the other metrics place a different emphasis on richness and evenness.
We also tested for statistical significance using paired Wilcoxon signed-rank test and Kruskal Wallis test.
After running estimate_richness, we tested for significance.
We found that water stress did not directly affect bacterial or fungal alpha diversity.
In bacteria, Shannon and Simpson indices decreased between control and drought in five of the six sites, but this was not significant (p > 0.05). However there were site-specific effects with significant differences between sites in the drought samples but not in the controls. Water deprivation can have multiple effects on diversity, some positive, some negative, depending on the literature. Positive effects have been attributed to adaptation and isolation reducing interspecies competition, while negative effects have been attributed to reduced nutrient availability, increase in osmotic stress and restricted mobility.
For fungi, differences in alpha diversity between control and drought were inconsistent with no significant differences between sites either.
This supports suggestions in the literature that fungi are less sensitive and more resistant to water stress than bacteria.
Post-hoc pairwise testing was carried out for the site and group variables using Dunn test as it is the recommended test to be run after Kruskal-Wallis:
As there were significant differences between sites in the drought samples, we carried out post-hoc testing to see where these differences occurred. This was only for the Shannon and Simpson metrics
Distances between samples were also calculated using ordinate method. The distance metrics selected were weighted and unweighted UniFrac for bacteria and Bray-Curtis and Jaccard for fungi.
Beta diversity analysis showed a clear separation between control and shelter samples for the unweighted UniFrac metric, which is sensitive to rarer taxa. This suggests that the effects of drought on rarer taxa was driving differences in community composition.
Additionally, control and shelter samples from the same site clustered together in most cases. However, this clustering was not as compact when compared to weighted UniFrac (see further down). This also highlights the impact of rarer taxa on community composition within site in response to water stress.
For weighted UniFrac, the PCoA plot showed no treatment-based sample clustering, while the effect of treatment was not statistically significant, indicating that short-term water stress had a stronger impact on rare bacterial taxa.
The pairing of samples from the same site was also more compact for weighted UniFrac when compared to unweighted.
We also carried out analysis using Bray-Curtis, but as these showed similar results to weighted UniFrac, they were not included in the overall analysis.
For fungi, differences in community composition related to treatment were not evident on the PCoA plots. There was no noticeable clustering of samples related to treatment for either Bray-Curtis or Jaccard along the top three PCoA plot axes. Shelter and control samples from the same site were paired on the plots although the closeness of samples was site-dependent.
Differences in community composition between samples based on treatment, site, and treatment-site were tested for statistical significance using PERMANOVA via vegan’s adonis2 function for all calculated distance metrics.
Statistical testing showed significant differences in bacterial communities between treatments for unweighted UniFrac. This indicated that water stress affected bacterial communities for rarer or less abundant taxa, supporting the PCoA plot for unweighted UniFrac
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = unifrac_dist.bac ~ Treatment + Site, data = metadata.bac, by = "terms")
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.05273 0.11294 1.6657 0.013 *
## Site 5 0.25586 0.54803 1.6165 0.001 ***
## Residual 5 0.15828 0.33902
## Total 11 0.46688 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = unifrac_dist.bac ~ Treatment + Site, data = metadata.bac, by = NULL)
## Df SumOfSqs R2 F Pr(>F)
## Model 6 0.30860 0.66098 1.6247 0.001 ***
## Residual 5 0.15828 0.33902
## Total 11 0.46688 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
However, once abundance was taken into account through weighted UniFrac, the differences in community composition between treatments were not significant.
There is support in the literature than repeated drying-rewetting events can alter taxonomic profiles of microbial communities through effects on rarer taxa. This has been attributed to rarer taxa having more specific resource requirements or just being more sensitive to environmental disturbance in general.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = wunifrac_dist.bac ~ Treatment + Site, data = metadata.bac, by = "terms")
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.003162 0.03713 0.8583 0.534
## Site 5 0.063599 0.74661 3.4523 0.002 **
## Residual 5 0.018422 0.21626
## Total 11 0.085184 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = wunifrac_dist.bac ~ Treatment + Site, data = metadata.bac, by = NULL)
## Df SumOfSqs R2 F Pr(>F)
## Model 6 0.066762 0.78374 3.02 0.001 ***
## Residual 5 0.018422 0.21626
## Total 11 0.085184 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = bray_dist.bac ~ Treatment + Site, data = metadata.bac, by = "terms")
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.009715 0.03179 0.7017 0.670
## Site 5 0.226687 0.74173 3.2749 0.001 ***
## Residual 5 0.069219 0.22649
## Total 11 0.305621 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = bray_dist.bac ~ Treatment + Site, data = metadata.bac, by = NULL)
## Df SumOfSqs R2 F Pr(>F)
## Model 6 0.236401 0.77351 2.846 0.001 ***
## Residual 5 0.069219 0.22649
## Total 11 0.305621 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
For fungal communities, differences in community composition were not significant between treatments for either Bray-Curtis or Jaccard. However, differences between sites, and between samples were significant. This suggests fungal communities varied greatly between sites while the impact of short-term water stress was modest and site dependent.
Factors which can affect fungal community composition include land management, soil pH and location. For example fungal communities can be negatively affected by intensive land management due to fertiliser application.
Comparing bacterial and fungal communities indicates a weaker effect of water stress on fungal taxonomic profiles. This reinforces the suggestion that fungi are less responsive to drought overall.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = bray_dist.fun ~ Treatment + Site, data = metadata.fun, by = "terms")
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.04056 0.05217 1.0789 0.423
## Site 4 0.58654 0.75441 3.9005 0.001 ***
## Residual 4 0.15038 0.19342
## Total 9 0.77748 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = bray_dist.fun ~ Treatment + Site, data = metadata.fun, by = NULL)
## Df SumOfSqs R2 F Pr(>F)
## Model 5 0.62710 0.80658 3.3361 0.004 **
## Residual 4 0.15038 0.19342
## Total 9 0.77748 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = jacc_dist.fun ~ Treatment + Site, data = metadata.fun, by = "terms")
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.08834 0.06436 1.0251 0.431
## Site 4 0.93962 0.68451 2.7256 0.002 **
## Residual 4 0.34473 0.25114
## Total 9 1.37270 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = jacc_dist.fun ~ Treatment + Site, data = metadata.fun, by = NULL)
## Df SumOfSqs R2 F Pr(>F)
## Model 5 1.02796 0.74886 2.3855 0.002 **
## Residual 4 0.34473 0.25114
## Total 9 1.37270 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1